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Abstract 


The deterministic transport code HZETRN was developed for research scientists and 
design engineers studying the effects of space radiation on astronauts and 
instrumentation protected by various shielding materials and structures. In this work, 
several aspects of code verification are examined. First, a detailed derivation of the light 
particle (A<4) and heavy ion (A>4) numerical marching algorithms used in HZETRN is 
given. References are given for components of the derivation that already exist in the 
literature and discussions are given for details that have been neglected in the past. The 
present paper provides a complete description of the numerical methods currently used in 
the code and is identified as a key component of the verification process. Next, a new 
numerical method for light particle transport is presented, and improvements to the 
heavy ion transport algorithm are discussed. A summary of various coding errors 
discovered while implementing the new method is also given, and the impact of these 
errors on previously predicted exposure quantities is shown. Finally, a coupled 
convergence study is conducted by refining the discretization parameters (step-size and 
energy grid-size). From this study, it is clearly shown that past efforts in quantifying the 
numerical error in HZETRN were hindered by single precision calculations and 
computational resources. It is also determined that almost all of the discretization error 
in HZETRN is caused by the use of discretization parameters that violate a numerical 
convergence criterion related to charged target fragments below 50 AMeV. Total 
discretization errors are given for the old and new algorithms to 100 g/cm 2 in aluminum 
and water, and the improved accuracy of the new numerical methods is demonstrated. 

Run time comparisons are given between the old and new algorithms for three 
applications in which HZETRN is commonly used. The new algorithms are found to be 
almost 100 times faster for solar particle event simulations and almost 10 times faster for 
galactic cosmic ray simulations. The improved algorithms will be implemented in all 
future versions of HZETRN. 

1. Introduction 

As human exploration beyond Earth’s orbit into radiation environments where measured data are 
sparse and testing is difficult, models will be heavily relied upon to make decisions regarding vehicle 
design and mission planning. This reliance on model results requires a systematic effort of verification, 
validation, and uncertainty quantification. Verification is the process of determining that a computational 
model accurately represents the underlying mathematical model and its solution; validation is the process 
of determining if the underlying mathematical model accurately represents physical reality, and uncertainty 
quantification is the process of identifying all relevant sources of uncertainties and quantifying their impact 
on the inputs and outputs of the model [NASA 2008]. This paper addresses verification of the radiation 
transport code HZETRN (High charge (Z) and Energy TRaNsport) [Wilson et al. 1991, 1995, 2006; Slaba 
et al. 2010] with a focus on documentation, improving efficiency and stability, and quantifying 
discretization error through convergence testing. 

Documentation is a critical component of verification [Roy 2005] and has been emphasized in the 
NASA standard for modeling and simulation tools [NASA 2008]. While there have been many papers 
published that describe the transport model and physical parameters [Wilson et al. 1975, 1977, 1986, 1991, 
2005, 2006] as well as the marching algorithms and numerical methods [Wilson et al. 1989, 1991, 2005, 
2006; Lamkin et al. 1992, 1994; Shinn et al. 1991] used in HZETRN, various gaps still exist in the 
documentation. For example, many of the mathematical details associated with inverting the Boltzmann 
integro-differential equation into a Volterra type integral equation have been neglected in the literature. 
Since the inversion depends on scaling conventions and physical arguments that are specific to the one- 
dimensional space radiation transport problem, these mathematical details are a necessary component of 
code documentation. In this paper, we first give a detailed review and derivation of the existing light 
particle (n, p, 2 H, 2 H, 3 He, 4 He) and heavy ion (A > 4) marching algorithms used in HZETRN. This will 
include all of the details associated with recasting the integro-differential equations as Volterra type 
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integral equations and all of the physical and numerical approximations. References are given for 
components of the derivation that already exist in the literature, and discussions are given for details that 
have been neglected in the past. As a result of the review, a numerical convergence criterion is identified 
that, to the authors' knowledge, has yet to be documented or examined. It is shown later in the paper that 
the discretization parameters (spatial step-size and energy grid-size) commonly used in HZETRN violate 
this criterion and cause a systematic under-prediction of light charged target fragments below 50 AMeV. 

This detailed derivation and review of the numerical methods also resulted in the development of a 
new light particle marching algorithm that is almost 100 times faster than its predecessor for solar particle 
event (SPE) simulations. Though computational efficiency has long been a core feature of HZETRN, there 
were certain applications for which the existing algorithms resulted in long run-times. For example, 
consider the interpolation or ray-by-ray methods used to compute mass averaged quantities in human 
phantoms exposed to space radiation. Interpolation methods are quite fast once the interpolation database 
has been generated, but it takes the current code over seven hours to generate a detailed database on a 
single processor. Similarly, the ray-by-ray method can take over 20 hours to compute the mass averaged 
particle fluence spectra at a single point in the body on a 192 core cluster. Approximately 1000 body points 
have been identified by Slaba et al. [2009] as sufficient for computing whole body effective dose in human 
phantoms exposed to SPE and galactic cosmic rays (GCR); this would indicate a significant computational 
cost. To help reduce these run-times and increase code efficiency, we present a new numerical method for 
the light particle marching algorithm that reduces the required number of interpolations and removes the 
need for integral fluence to be calculated at each step. The new method is shown to be almost 100 times 
faster for solar particle event (SPE) simulations and almost 10 times faster for galactic cosmic ray (GCR) 
simulations. The accuracy of the new methods is discussed later in the report. 

Controlling round-off error and identifying coding errors is another important point of code 
verification. Though previous convergence studies [Shinn et al. 1991; Slaba et al. 2010] and benchmark 
comparisons [Wilson et al. 1988, 1991, 2005, 2006; Heinbockel et al. 2009a, 2009b] would indicate that 
round-off error is of little concern, such comparisons were generally made at moderate shielding depths 
where round-off errors are assumed to be small. However, as HZETRN is increasingly used in atmospheric, 
lunar, or Martian surface applications with large material thicknesses (>50 g/cm 2 ), round-off error could be 
a major concern and needs to be investigated. Therefore, we also discuss the impact of round-off error and 
coding errors in the existing marching algorithms in HZETRN. Selected light particle cross sections are 
calculated in single and double precision, and the impact of round-off error in the single precision 
calculations is shown to be large in certain cases. In the process of modifying HZETRN to enable double 
precision calculations, various coding errors were also discovered that have non-negligible effects on 
particle fluence spectra, dose, and dose equivalent. The coding errors are discussed and resolved. The 
interpolation routine [Wilson et al. 1995] frequently used in the transport algorithms is also examined, and 
a new routine is developed that is faster, has improved extrapolation procedures, and has the capability of 
interpolating around certain discontinuities. The improved code stability attained by using double precision 
calculations and removing the various coding errors is clearly shown. 

In many computational models or algorithms, continuous variables are discretized to reduce a 
differential equation into an algebraic expression that is evaluated numerically. The algorithm is said to 
converge if the numerical solutions reach an asymptotic limit as the discretization parameters approach 
zero. In order to show convergence and quantify discretization error, the discretization parameters are 
refined several times and the differences between the various solutions are compared. Such studies are 
often referred to as convergence tests. As part of a larger verification and validation effort, configuration 
controlled convergence tests are created which can be rerun when significant changes are made to the 
codes. The ability to rerun such tests will help prevent the introduction of errors into the code as 
modifications are made in the future. In HZETRN, the spatial variable x and energy variable E are 
discretized. Two convergence tests have previously been published [Shinn et al. 1991; Slaba et al. 2010]; 
however, those tests were primarily focused on verifying code stability and were limited by computational 
resources. In the first analysis, Shinn et al. [1991] conducted a coupled convergence test in both space and 
energy. However, only two step-sizes (h = 0.5 g/cm 2 and h = 1.0 g/cm 2 ) and two energy grids ( N = 30 and 
N = 60, where N is the number of grid points) were considered; only nucleons were transported, and all 
calculations were in single precision. No attempt was made to quantify discretization error, and numerical 
convergence was not clearly demonstrated. In the more recent analysis, it was determined that at least 100 
energy grid points are needed to control energy discretization error [Slaba et al. 2010]. It will be shown 
here that even step-sizes of h = 0.01 g/cm 2 can result in moderate errors for low energy target fragments 
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with small residual ranges. Neither the codes nor the results were configuration managed in these previous 
studies. Thus, a more detailed convergence analysis is necessary. We conduct a convergence study in both 
step-size and energy grid-size. The spatial discretization parameter, h , is reduced by factors of 2 from its 
common value of 0.5 g/cm 2 down to 2' 11 g/cm 2 . Similarly, the number of energy grid points is increased by 
factors of ~1.5 from 100 up to 753. Particle fluence, dose, and dose equivalent values are then computed at 
various depths in aluminum and water slabs exposed to SPE and GCR environments using the old and new 
transport algorithms and each of the discretization parameters. The resulting data are used to show that the 
new algorithms reach an asymptotic solution as the discretization parameters are refined. The improved 
accuracy of the new methods is also clearly demonstrated. Discretization errors are also given for the 
discretization parameters commonly used in HZETRN (h = 0.5 g/cm 2 with 100 energy grid points). These 
errors are expressed as percent differences from the converged numerical solutions obtained with the finest 
discretization parameters. 

2. Transport Equations 

The Boltzmann transport equation with the continuous slowing down and straight ahead 
approximations is given as [Wilson et al. 1991] 

B[^{x,E)} = ?/; a jk (E,E')<l> k (x,E')dE', (1) 


with the linear differential operator 


B[^(x,E)\ 


d_ 

dx 


3 


(j>j{x,E), 


( 2 ) 


and boundary condition 


^•(0 ,E) = f 3 (E). 


(3) 


In equations (l)-(3), (j>j(x,E) is the fluence of type j particles at depth x with kinetic energy E , Aj is the 

atomic mass number of a type j particle, Sj(E) is the stopping power of a type j ion with kinetic energy E, 
Oj{E) is the total macroscopic cross section for a type j particle with kinetic energy E , and gj^E^E) is the 
macroscopic production cross section for interactions in which a type k particle with kinetic energy E 
produce a type j particle with kinetic energy E. The summation limits in equation (1) will be discussed 
shortly. The boundary condition spectrum, f. (E ) , is considered to be a known function over a broad 
energy spectrum. 

Consider the continuous slowing down operator 


— — —S-(E) 
A. dE 3 ’ 


<l>i(x,E), 


(4) 


which represents the rate at which charged particles lose energy as they interact with the electron clouds in 
the target media. Although atomic interactions cause charged particles to lose energy in discrete increments 
as they pass through a material, there are a sufficient number of these interactions in a unit path length to 
justify a continuous approximation [Wilson et al. 1991]. It is advantageous to approximate this term by 
considering the relation from Bethe stopping power theory [Bethe et al. 1930] 


v r . , 

3 3 9 


(5) 
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where v. = 2?. / A. with Zj being the atomic charge of a type j particle, and the range of a type j ion, r 7 -, is 
given by the range-energy relationship 


A.f'JS-. 


From equation (5) and (6), it is clear that the proton range and stopping power can be approximated as 


r 


VT- 3 
3 3 5 


(7) 


( 8 ) 

3 

The approximation in equation (7) is depicted in Figure 1 and is less accurate at low energy (<10 AMeV) or 
low residual range. 

Equation (8) allows the transport equation to be written as 


with 


B[^(x,E)} = £ 
k 


a jk (E) E ')<p k (x, E')dE' , 


(9) 


B[^(x,E)} 


d 

V 

dx 3 






( 10 ) 



Figure 1. Comparison of scaled proton range, r/vj, and 4 IIe range ,r } , in aluminum. 
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The motivation, justification, and physical implications associated with approximating each ion stopping 
power with a scaled proton stopping power has been discussed in detail elsewhere [Wilson et al. 2006]. 
Ultimately, the approximation allows significant numerical simplifications, as will be shown later. 

There are two paths taken (for light particles and heavy ions) in developing numerical procedures 
for equation (9). For heavy ions, it is noted that projectile fragments have energy and direction very near 
that of the projectile, while target fragments are produced nearly isotropically with low energy and travel 
only a short distance before being absorbed [Wilson et al., 1977, 1986, 1991, 1995, 2006]. The approximate 
decoupling of target and projectile fragments is discussed in detail by Wilson et al. [1977, 1991, 1995] and 
suggests that target fragments can be neglected in the heavy ion transport procedure (their contribution to 
dose is approximately accounted for after the transport procedure). The equal velocity assumption for 
heavy ions can be expressed in the production cross section as [Shinn et al. 1992] 

<r jk (E,E') = a jk (E')6(E-E'), (11) 


where ojilE) is the production cross section for interactions in which a type k particle with kinetic energy E 
produce a type j particle. In this case, equation (9) becomes 

= ( 12 ) 

k 

The absence of target fragments in the heavy ion transport procedure allows one to take the summation in 
equation (9) over all projectiles with mass greater than that of the fragment. If all the transported particles 
are ordered according to mass, then equation (9) can be succinctly written as 

B[<t> j (x,E)]='£ / a jk (E)<t> k (x,E), (13) 

k>j 

which is the transport equation found in Wilson et al. [1986, 1991, 1995, 2006] and will be referred to as 
the heavy ion transport equation. The upper summation limit in equation (13) can vary, but it is common to 
use no fewer than 59 ions. See Cucinotta et al. [2006] for a discussion of isotope selection. 

Alternatively, for light particles (n, p, 2 H, 2 H, 3 He, 4 He) both projectile and target fragments are 
included in the transport procedure. The broad energy distribution in collision events also indicates that the 
equal velocity assumption in equation (11) cannot be used. In this case, no simplifications to equation (9) 
are used, and the summation is taken over all light particles. Hereafter, equation (9) will be referred to as 
the light particle transport equation which includes both neutrons and light ions. It should be noted that for 
SPE environments with a negligible heavy ion component, only the solution to the light particle transport 
equation is required. For GCR environments, there are both heavy ion and light ion components, and 
solutions to the light particle and heavy ion transport equations must be evaluated simultaneously. The 
coupling of these two equations will be discussed later. 

3. Light Particle Marching Equation 

The following formulation of the light particle marching equation was taken from Wilson et al. 
[1975, 1989, 1991, 2006], Shinn et al. [1991], Lamkin et al. [1992, 1994], Cucinotta [1993], and Slaba et 
al. [2010]. Further references in this section will be given only for mathematical techniques and physical 
arguments. 

To develop a numerical marching procedure for the light particle transport equation, it is necessary 
to invert the differential operator on the left hand side of equation (9) by using the method of characteristics 
[Haberman 1989]. To do this, define the scaled quantities 

= UjS(E)<j>j(x, E) (14) 


and 
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s jk (r,r') = S(E)a jk (E,E'), (15) 

where the proton range, r, is defined in terms of the proton stopping power, S(E), by the range-energy 
relationship 


dE ' 

0 S(E') 9 


and the modified scaling parameter is given by 


1 , j = n 
Vj , j * n' 


(16) 


(17) 


The quantity v . is used in fluence scaling to avoid a trivial solution for the neutrons. Notice that we have 

used the scaled proton range in equation (14) to allow the function tjj. {x, r) to be defined over the common 
proton range, r, for each particle type j. 

In terms of the scaled quantities introduced in equations (14) and (15), the integrals on the right 
hand side of the light particle transport equation become 


e/; <T jk (E,E')<t> k (x,E')dE' = e/: 


1 8(B) 


S j k ( r > r ') 


lw 


ip k {x,r') 


dE' 


E f r °° s jk (r,r')i> k (x,r') dr' , 


(18) 


where the differential has been transformed using 


dr' 


S(E ') 


-dE' 


(19) 


Now, note that 


fol>i _ foPj dE 
dr ~ dE dr 


u.S(E)-^[S(E)<t,.(x,E)], 


( 20 ) 


so that the differential operator on the left hand side of the light particle transport equation can be expressed 
in terms of the scaled quantities as 


d 

dx Vj 


dE 


S(E) + *.(E) 


4>Ax,E) 


a+L 

dx 



s (E)<t>j] + < JjiEty 


1 dtpj Vj dtp j 1 

VjS{E) dx ~ v 3 S{E) dr + v f S(E) 


d 

dx 


Tr + 



( 21 ) 


where the functional dependencies of the scaled and unsealed fluence have been suppressed in certain 
places for brevity. Equations (18) and (21) are now combined to give the light particle transport equation in 
terms of the scaled quantities as 
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( 22 ) 


--V. — + oAE) 
dx 3 dr 3 


1 . ^ 1 poo 

ibAx.r) = > I s- t (r,r'm t (x.r')dr 

v-S{E) 3 . T7. ft( E}Jr > kK ’ >YkK ’ 3 


* k S{E) 


which simplifies to give the final scaled form 

d 


d_ 

dx 


--V. — + aAr) 
Q 3 dr 3 


k ‘'k 


^ ^ V ■ poo 

& r) = ~ J r s jk ( r , r ')tp k (x, r')dr' , 


(23) 


The term a ■ (E) has been replaced with a- (r) , since for a given value of r, equation (16) can be inverted 
to determine E. 

Equation (23) is inverted (see Appendix A) using the method of characteristic to obtain the 
Volterra type integral equation 


ipfar) = e + v.x) 


+V'z L f f e ^ r ^ s j k { r + VjX\r')il; k (x - x\r ')dr'dx' , 

V , J 0 Jr+vp' J J 


(24) 


with the integrating factor 


0j(r,x) = f Q a-(r + v.t)dt . (25) 

Equation (24) can be easily written as a marching procedure in terms of the step-size, h , as 

ipjix + h,r) = e-fy^ip^x.r + vji) 

, — f'/l poo of (26) 

+y j , e r ’ 2 s jk( r + } )^ k {x + h — x\r')dr y dx y . 

& T+VjH 

We now seek to approximate the differential fluence in the integrand of equation (26). Define 

Q jk (x,r + v.x ') = J + e -^( r ^’) 5 ^( r + + h - x\r')dr 1 , (27) 

then 

<3# (z, r + ') = Q jk (x, r) + 0(x ’) . (28) 

The 0{x k ) notation indicates that the remaining terms in the polynomial expansion are of degree k or 
greater. The exponential attenuation term can also be expanded as 

e -i 9 fa') =i + 0 {x'), (29) 

where the expansions e~& = 1 — f3 + 0(/3 2 ) and (3j(r y x } ) « x'<Tj(r) have been used. Substitution of 
equations (27)-(29) into equation (26) produces 
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( 30 ) 


ipjix + h,r) = e + Vjh) + zr Q jk (x, r) + 0{x')dx 

k V k 


or 


ipjix + h,r) = e ip far + v.K) + 0(h ) . 

Equation (31), without the 0(h) term is the exact solution to the homogeneous equation 


(31) 


9 9 _L 

V 4 h (TAT ) 

3 dr 3 


4>jfar) = 0 . 


(32) 


The homogeneous equation neglects secondary particle production through nuclear interactions and 
accounts only for the slowing down of particles due to atomic interactions and the loss of particles due to 
nuclear absorption. If the step-size is taken to be sufficiently small such that 


h < 



(33) 


(i.e. much less than the nuclear mean free path), then the local truncation error will be negligible as the 
particles will not have travelled far enough to suffer a nuclear collision [Wilson et al. 1991]. Equation (33) 
is the first convergence criterion and has been well documented in the literature [Wilson et al. 1975, 1989, 
1991, 2006]. It is worth noting that nuclear mean free path lengths are on the order of many g/cm 2 , while 
the step-sizes usually taken in HZETRN are less than 1 g/cm 2 . 

Equation (31) can be used to approximate the integrand in equation (26) by noting that [Wilson et 
al. 1975, 1989, 1991,2006] 

tp k (x + h — x\r') = e~^ r ' ,h ~ x ^ip k (x : r } + v k (h — x')) + 0(h) . (34) 


Substitution of equation (34) into equation (26) yields 


tl>j(x + h,r) = e ^ (r ’^(a,r + vfr) 
v ■ 


^ V k '*Q J r+1/jX' 

+o(h 2 ) , 


)-/3 k (r\h,-x') 


s jk (r + is j x\r')ip k (x,r'+ v k (h — x'))dr'dx'( 35) 


where the 0(h 2 ) terms are small. This is the light particle marching equation given by Wilson et al. [2006] 
for which numerical procedures will be developed and studied later. Hereafter, the 0(h 2 ) will be assumed 
in the marching equations and not written. 


4. Heavy Ion Marching Equation 

The following derivation of the heavy ion marching equation was taken from Wilson et al. [1977, 
1986, 1991, 1995, 2006] and Shinn et al. [1992]. Further references will only be given for mathematical 
techniques. Any variables and symbols used in this section that were defined in the previous section hold 
their respective definitions. 

The heavy ion transport equation was given previously in equation (13) as 


8 



( 36 ) 


B[(/>.(x,E)] = ^a jk { E )<i>k{x,E). 

k>j 


Define the scaled heavy ion fluence 


= v-S(E)4>.(x,E), (37) 

and replace a jk ( E ) with a jk (r) since, for a given r, equation (16) can be inverted to find E. Notice that the 
scaling in equation (37) is identical to the scaling used for the light particles in equation (14) except that we 
have used v. instead of v. . This minor change reflects the fact that we are now only transporting heavy 

charged particles for which v j >0 and trivial scaling is not possible. 

Following the procedure outlined in the previous section, equation (36) is written in terms of the 
scaled fluence as 


a d 

Va b cr.fr ) 

dx 3 dr 3 


v . 

ipjM = 22 —<r jk (r)i(> k (x,r). 


k>j ^ k 


(38) 


Equation (38) is inverted (see Appendix B) using the method of characteristic to obtain the Volterra type 
integral equation 


ipj(x,r) = e ,r + vp) 

+y^J 0 (r + v.x 1 ) Ip k (a: - x\r + v^dx' , 

k>j V k 

with the exponential term 

/3j(r,x) = J'J (Tj(r + Vjt)dt . 

Equation (39) can be easily written as a marching procedure in terms of the step-size, h , as 

+ h y r) = e~^ r ^fpj(x 9 r + isfi) 

. V. nh 


(39) 


(40) 


(41) 


(42) 


J e ^ r ' x ^(7j k (r + v.x y ^ r ip k i^x + h — x\r + v^x'^dx' . 

k>j 13 k 

As before, use the 0(h) homogenous solution [Wilson et al. 2006] 

ip k (x + h — x\r + VjX ’) = e~^ r+I/jX ’ h ~ x ^(a;,r + v.x '+ v k (h — x ’)) + 0(h ) , 

to obtain 

ipj(x + h,r) = e~ /3 > ( - r ’ h ' ) 'ip j (x,r + vh) 

+y^f G e _ ^ r,x ')-Pk( T + v f c < h ~ x _)_ v .x'^tp k (x,r + v.x'+ v k (h — x'))dx' (43) 


k>j k 

+0(h 2 ) . 
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From here, approximate the integral over x ' from 0 to h by expanding the production cross 
section and particle fluence in a Taylor series [Wade 2000] as 

a jk ( r + v i x ') = a jk ( r + v j h / 2 ) + °( x ')> ( 44 ) 

^k(x,r + v-x '+ v k (h - x'))= if> k (x,r + (v j + v k )h / 2) + 0{x') , (45) 

and obtain 


i>j(x + h,r) = e ^ r,h ^7p.(x,r + v.h) 

+ J2~ a jk ( r + v j h / 2 )A( x ’ r + i v i + v k) h / 2) \j\^ {r ’ x ' ) ~ 0k(r+v i x '’ h ~ x ' ) dx ' (46) 

k> j V k 

+ 0(/i 2 ). 

The remaining integral is similarly handled by considering the following Taylor expansions for the 
exponential terms 


Pj(r,x') = JJ a-{r + vfidt 

~ fj a j( r + v 3 h / 2 ) dt 

= a-(r + v-h / 2)z' + o((a:') 2 ), 

P k (r + VjX',h -x') = f* X a k (r + v.x'+ v k t)dt 

~ fo X a k( r + u j x '+ v k h / 2 ) dt 

= a-(r + v-x '+ v k h / 2) (h - x ') 

~ <T ji r + i v j + v k ) h / 2 )( h “ *') + ^((a:') 2 ) • 


(47) 


(48) 


The integral is now evaluated as 


C h e -Pj{r,x')-P k (r+u j x\h-x') d x y ^ C h e -<T j (r+i/ j h/2)x'-<r k (r+(u. 

Jo ~ Jo 


j +u k )h,/2)(h-x') 


dx 1 


a k( r + ( v j + v k) h / 2 ) - a j( r + v j h / 2 ) 


+ 0(h 2 ) . 


(49) 


The final marching procedure is given by 

ipj(x + h,r) = e^^ip-far + v.h) 

+ ( r + v j h / 2 )M x ’ r + ( y, + v k) h / 2) 

k>j V k 

-Vjir+Vjh/tyh _ -cr k {r+{v j +v k )h/2)h 

X + 0(h 2 ) , 

a k( r + ( v j + v k) h / 2 ) - <Tj(r + Vjh / 2) 


(50) 
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where the 0(h 2 ) terms are small. This is the heavy ion marching equation given by Wilson et al. [2006] for 
which numerical procedures will be developed and studied later. Hereafter, the 0(h 2 ) will be assumed in 
the marching equations and not written. 

5. Existing Light Particle Marching Equation Numerical Methods 

In this section, the numerical methods developed to evaluate equation (35) are reviewed and the 
error analysis found in the literature for those methods is summarized. The following notation and 
terminology will be used extensively throughout the remainder of the report. An energy grid (E-grid) refers 
to a discrete set of energy values distributed in some manner between a minimum energy value, , and 

a maximum energy value, E max . The z' th component of the E-grid is denoted as E { . A range grid (r-grid) 
refers to a discrete set of proton range values distributed in some manner between a minimum range value, 
Tallin , and a maximum range value, r max . The z ,th component of the r-grid is denoted as r. . The number of 

grid points in a grid, or grid-size, is denoted by N. It will be assumed that all grid indexing is from z = 1 to 
i = N . Equal-log spacing is also used extensively; the f 1 component of an equal-log spaced E'-grid is 
evaluated as 


= 20 lQg ( jB ndn)+ A B(*- 1 ) 9 

where 


lQ g (£ m J ~ jggCfinm) 

N -1 


(51) 


(52) 


and the logarithms are base 10. Note that an energy grid can be converted to a range grid, and vice-versa, 
using equation (16). 

In Section 3, we derived the light particle marching equation (35) as 


ipjix + A,r) = e ^' (r,A Vj(z,r + v jh) 

+ ^zrf Q f + ,h ~ x ^s jk (r + v .x\r')ip k {x,r '+ v k (h - x'))dr'dx\ 


(53) 


It is clear from equation (53) that the quantity tjj- (x, r + v.h) will likely have to be evaluated for range 

values not on the prescribed r-grid. Wilson et al. [1991] and Lamkin et al. [1992] have shown that a log-log 
cubic Lagrange interpolating polynomial provides sufficient accuracy with a minimal number of grid 
points. The terminology “log-log” refers to taking the logarithm of both the dependent and independent 
variables prior to interpolating. The interpolated value is then exponentiated to adjust for the initial 
logarithm. It has been shown that higher order Lagrange polynomials offer little more in accuracy, and 
higher order splines can introduce uncontrollable oscillatory behavior that is problematic in iterative 
marching procedures. 

First simplify the integral over x 1 from 0 to h by using a modified single point midpoint rule for 
integrable functions a( x) 9 b(x), and c(x) of the form 


X h 

a(x)b(x)c(x)dx 


<h/ 2) /: b(x)dx 


c(h / 2) . 


(54) 


This approximation will be accurate if the functions a(x) and c(x) are nearly constant over the interval 
[0 ,h \ . When applied to equation (53), this approximation yields [Wilson et al. 1991, 2006] 
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ipj{x + h,r) = e + v.h) 


+ 


V a poo 


E I-! 

* p t Jr 


-^h/2)-^(r\h/2) \ rh 


r+vJi/2 


f 0 s jk( r + vp',r')dx' A( x > r '+ v k h / 2 ) dr ' 


( 55 ) 


It is clear that the modified midpoint rule has decoupled the source integrals in equation (53); however, the 
approximation will only be accurate if the step-size is sufficiently small to satisfy 


e(r + v.h / 2) w e(r) , 


(56) 


or 


h < r / v. « r.. (57) 

Equation (57) is the second convergence criterion. Recall that the first criterion was identified 
previously in equation (33) and required the step-size to be smaller than the nuclear mean free path length. 
To the authors' knowledge, this second criterion has never been explicitly stated or addressed in the 
literature. Note that the range of a 100 MeV proton in aluminum is -10 g/cm 2 , so a step-size of 0.5 g/cm 2 
will be sufficient in this approximation. However, a 10 MeV proton in aluminum will travel -0.17 g/cm 2 
before coming to rest, and a 1 MeV proton in aluminum will travel -0.004 g/cm 2 before coming to rest. 
This indicates that step-sizes of 0.5 g/cm 2 and 1.0 g/cm 2 studied previously [Shinn et al. 1991] do not 
accurately transport low energy target fragments. It will be shown later that even step-sizes near 0.01 g/cm 2 
will systematically under-predict the source integral by neglecting the particle production from projectiles 
with low ranges near that of fragments. This approximation will be tested in detail when the coupled energy 
grid and step-size convergence study is conducted. 

Continue simplifying the light particle transport equation by noticing that for charged particle 
fragments 



( \ /»r+i/ 'Jh 

v-x r ') dx 1 = — I s ik (u, r ') du 

d i / , J r d 

J 

1 [ nr+v-h nr 

= ~\J 0 a jk( u ’ r ') du - J 0 *#(«>*•')<* 

= v:\C w,l ‘ >,T * iE ’ E ' )dE ~ 


F£(r,r']h) . 


(58) 


The term e(r) is the energy associated with the proton range r, and for neutron fragments, denoted by j = 
n 9 



v n x\r')dx 1 



= K*( r > r ') 

= F&(r,r'-,h) . 


(59) 


The light particle transport equation is now written as 


ipjix + h,r) = e 0 >^'ip j (x,r + vti) 


f°° h/2 e ^ r,h/2) 0k{r ' ,h/2) F£(r,r']h)ip k (x,r'+ v k h / 2) dr' 

k V k Jr+ Vj h/2 


(60) 
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The exponential attenuation terms are simplified using [Shinn et al. 1991; Larrikin et al. 1992] 


and 


tyrfi) = e -fo a > (r+V > t)dt ~ 


nhj l 

e -0 } (r,h/2) = e ~J 0 <Tj(r+Vjt)dt _ e -<7 j (r)ft/2 


ntlfl 

e -p t (r\h/ 2) = e ~J 0 ^(r'+v&dt _ e -a k (r')hl2 


( 61 ) 

(62) 


(63) 


Equations (61)-(63) are exact for neutrons and have been shown [Shinn et al. 1991; Lamkin et al. 1992] to 
produce negligible errors for charged particles with step-sizes up to h = 0.5 g/cm 2 . The light particle 
transport equation is now reduced to 


ipjix + A,r) = e <Tj ^ h tf) : j(x,r + v^h) 


I V 3 f 00 -<7j(r)h/2-a k (r')hl2 

* 


(r, r'-,h)ip k (x, r '+ v k h / 2 )dr' . 


(64) 


Now consider another one-point quasi-midpoint rule for integrable functions f(x) and g(x) 
integrated over a closed interval [a, b\ 


f f(x)g(x)dx & f(x) [ b g(x)dx , (65) 

a i/ a 

where x = (a + b) / 2 . This approximation will be accurate if f( x) is nearly constant over the interval 
[a, b\. The composite quasi-midpoint rule is obtained by extending the approximation over several sub- 
intervals within a given region 


n b 

/ f(x)g(x)dx 

a 


D /(*<), , 


( 66 ) 


where = a , x N = b , and a;. = (z. +1 + z.) / 2. This approximation is similar to the multi-group 

method commonly used in nuclear reactor transport theory and requires the sub-intervals [x v x i+1 \ to be 

sufficiently small such that f( x) is nearly constant over the sub-interval [Marchuk et al. 1986]. The 
composite quasi-midpoint rule in equation (66) can be applied to the integral in equation (64) by first 

making a simple change of variables (r" = r — v.h / 2) so that the light particle transport equation 

becomes 


ipjix + h,r) = e + v jti) 


+ 


Y]zr f e aj(r)A/2 +v 3 h l 2 ^ h l 2 F^ (r, r H + v.h / 2]h)'tp k (x,r ■■+ (v, + v k )h / 2) dr" 

U Vj. J r 3 3 J 


(67) 


Define the integral fluence as 
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( 68 ) 


poo 

&j(x,r) = Jr i/)j(x,r")dr" , 
so that for a given r h equation (67) is evaluated as 


ipfa + = e <Tj ^ h rp : j(x,r i + v-K) 

77 n- 1 


j \ 1 i ' j 


+E^ E (r„r„ + / 2;ft) 

fc m—i 

X l^j( x i r m + / 2) — ^j( x > r m+l + (^+^)ft/2)], 


(69) 


where r m = (r m+1 + r m ) / 2 . Slaba et al. [2010] have recently shown this approximation to be poor for 
neutron elastic interactions in which the energy loss between the pre-collision and post-collision neutron is 
very small. They have adjusted the neutron elastic component, F*f l , so that it is now evaluated as 


Fnn e %,r m ,h) = hSW^^rj), 


(70) 


where 


( a K r i> r m)) 


< r m+l) ~ £ ( r m) J£(r » ) 




(71) 


is the average value of the differential neutron elastic production cross section over the interval [r m , r m+ 1 ]. 

Note that another approximation has also been made, namely 

e -^(r n +^/2)fc/2 _ e -a k (r n )hl2 ' ^ 

Equations (68)-(72) define the light particle marching procedure used in HZETRN and tested in this report. 
A final point to make about the light particle marching procedure is related to manner in which it is coded 
(ie. the main algorithm used to evaluate the equation). Figure 2 gives a general process description of the 
light particle marching algorithm. The required inputs are the initial and final positions as well as the 
boundary flux/fluence at the initial position. The next step is to determine the maximum step-size (h max < 
0.5) and integer P such that P iterations of the light particle marching equation will propagate the fluences 
from the initial position to the final position. Once the constant step-size has been determined, the 
exponential attenuation term and production cross sections are computed for all values of r on the r-grid. 
These values are then used repeatedly within the main loop. Once the main loop is entered, the integral 
fluence must be computed, and the differential and integral fluences must be interpolated for each particle 
and each range value. 

The most computationally expensive portion of the calculation is the cross section calculation, 
which has been outlined in a bold dashed line in Figure 2. The integral fluence, interpolation, and source 
term calculations are not as computationally expensive, but are repeated several times within the loop and 
are therefore important. They have been outlined with a smaller dashed line. It should also be noted that in 
this algorithm, all of the computations are carried out in single precision. More will be said about this later. 
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6. New Light Particle Marching Equation Numerical Methods 


In this section, new numerical methods are presented to evaluate the light particle marching 
equation. The approximations used in the previous section to get from equation (53) to equation (67) are all 
used in this section; therefore, begin with equation (67) and make the substitution 

r* = r"+ (y. + v k )h / 2 , (73) 


so that the light particle transport equation becomes 


+ h,r) = e + vh) 



(r)h/2-a k (r'-v k h/2)h/2pA 

r jk 


(r,r — v k h / 2;h)xp k (x J r } )dr 1 , 


(74) 


where 



(75) 


Now, suppose that the fluence in the integrand, ip k (x, r ') , can be approximated by a linear 
combination of log-linear basis splines as 
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(76) 


N 

M X ’ r ') ~ Z B m ( r ')A (*> r m ) . 


m=l 


with the log-linear spline 




/ r m-l)/ ln ( r m / 

ln (r ro+ i / r )/Hr m+ 1 

0 



» r G [ r ra -l> r J 

» r e [ r m> r m+l] 
, otherwise 


(77) 


Substitute equation (77) into equation (74) and get 

^•(a: + ft,r) = e _<Tj(r)A ^.(a:,r + u.h) 

+ Z Z) +„ t )A /2 e~ a ’ (r)h/ 2 ~ at(r '~ Vth/ 2 )h/ 2 F^ (r, r / 2; ft)5« (r ') dr ' 1 . 


(78) 


ft 771=1 


It is clear from equation (78) that the integrand and the integral no longer depend on the fluence or the 
depth in the material. In fact, the integral can now be treated as a constant matrix of production coefficients 
depending only on the step-size and the r-grid. If we evaluate equation (78) at the z th r-grid value, then the 
production coefficients can be defined as 


.<*> s x 


j(*i)h/2-<r k (r'-v k h/2)h/2 


r i+( v j+ v k) h / 2 


F jk ( r i » r V k h / 2 ; h ) B m ( r ') dr ' 


so that the transport equation is now greatly simplified as 


(79) 


i’j (x + h, r.) = e °> {r ' )h ^ j (*, r. + v.h) + Z Z a im ( h )A (*> r m ) • ( 80 ) 

ft 771=1 

This approach casts the production integral as a linear combination of terms computed on the original r- 
grid requiring no further interpolation or integration. 

The production coefficients, a im , are further simplified by noting that the basis spline, B® ( r ') , is 
non-zero only in the region [r* m _ 1? r m+1 ] • Therefore, 


.« - x 


" b e ^j( r i) h / 2 -°k( r '- v k h / 2 ) h / 2 ] 


F jk( r i> r '~ v k h / %h)B®(r')dr' 


(81) 


where 


{ r « + ( v j + v k)h / 2 > r ro -i } , 

(82) 

= min { r N , r m+1 } , 

(83) 


and min and max refer to the minimum and maximum values, respectively. For neutron elastic interactions, 
the upper limit has been adjusted to 


r b = min K>w r a }> 


(84) 
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where r a is the proton range associated with the energy E a = £r(rj) / a T , efc) is the energy associated 
with the proton range r u and a T is the target dependent parameter given by 


Aj. -1 


Ap + 1 


(85) 


where A T is the mass of the target nucleus. To understand where the term r a was taken from, note that 

Wilson et al. [1991] have parameterized the neutron elastic production cross section to be non-zero in the 
energy region 


a T E' < E < E\ (86) 

where E is the energy of the pre-collision neutron and E is the energy of the post-collision neutron. This 
equation can also be expressed in terms of the pre-collision neutron energy as 

E < E' < E / a T . (87) 

For heavy targets, in which a T -+ 1, the non-zero region of the neutron elastic production cross section can 
become smaller than the spacing of the E-grid (or r-grid). Including the upper bound of the elastic cross 
section in the selection of the upper limit of integration will properly account for such occurrences. 

Similarly, for quasi-elastic interactions involving 4 He projectiles, the lower limit of integration has 
been adjusted to 


r a = maX { r i + ( V j + V k) h / 2 > r m-V r i ,3) + h / 2 } » ( 88 ) 

where is the proton range associated with the energy e(rj) + E B , and E B is the binding energy of the 

struck nucleus. To understand where the term + h / 2 was taken from, recall that for 4 He quasi-elastic 
interactions, secondary 4 He ions have energy E satisfying 

E<E'-E B9 (89) 

where E is the projectile energy and E B is the binding energy of the struck nucleus. In equation (81), the 
cumulative production cross section, , is evaluated at r y — v k h / 2. Since v k = 1 for 4 He, the 

minimum projectile energy must be e(r^ + h / 2) . 

In order to further simplify the integral in equation (81), the exponential terms in the integrand are 
approximated as 


e -v j (r i )h/2-(T k (r'-is k h/2)h/2 ^ e -ff i (r i )ft/2-o' Jb (r)ft/2 


(90) 


with r = (r 0 + r h ) / 2 . The approximation is valid since the exponential terms are slowly varying over the 
small interval (r a9 r b ). In this case, the production coefficients become 


,(*) = e -<T,(r<)ft/2-0 * (F)ft/2 f£ ( r { , r v k h / 2-,h)B®(r')dr' 


(91) 
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and are numerically evaluated using a 10-point Gauss-Legendre quadrature. For some reactions in which 
the cross sections are highly peaked or rapidly varying over the interval (, r a , r b ) 9 several subintervals are 
used to improve accuracy, and the quadrature is applied over each of the subintervals. For n producing n 
through elastic collisions, ten subintervals are used (100 quadrature points per integral); for n producing n 
through reactive interactions, five subintervals are used (50 quadrature points per integral); for 4 He 
producing 4 He through quasi-elastic interactions, five subintervals are used (50 quadrature points per 
integral). For all remaining reaction channels, two subintervals are used (20 quadrature points per integral). 

Even though the numerical integrations are carried out over small subintervals of the r-grid, the 
production coefficients are still computationally expensive to calculate for a single h value; therefore, they 
are approximated by using the log-linear interpolation 

a i m ( h ) * ex pj a !!(Vi) + [°£(Vn) - a S( ft p)]( ft(0 - h p ] ) / ( h iii - A ®)} = °L » ( 92 ) 

where = log(a im ) and ft® = log (ft) . This requires the coefficients a im (h) to be pre-computed for 
several values of h from 0.0 g/cm 2 to 0.5 g/cm 2 . The ft-grid has been taken to be 

h p = {0.0, 0.001, 0.005, 0.01, 0.05, 0.1, 0.25, 0.5} . (93) 

The final marching equation for light particles is now written as 


+ h,n) = e + v 3 h) + a imM x ’ r m) ■ ( 94 ) 

k 771=1 

After reviewing the older light particle marching algorithm depicted in Figure 2, it was determined that the 
cross sections used in the exponential attenuation term could also be pre-computed and passed to the 
algorithm as an argument. This modification has been made in the new algorithm. A general process 
description of the revised light particle marching algorithm is depicted in Figure 3. The required inputs are 
the initial and final positions as well as the boundary flux/fluence at the initial position, the attenuation 
cross sections, and the production coefficients as a function of h. The next step is to determine the 
maximum step-size (ft max < 0.5) and integer P such that P iterations of the light particle marching equation 
will propagate the fluences from the initial position to the final position. Once the constant step-size has 
been determined, the exponential attenuation term and the production coelficients are computed for all 
values of r on the r-grid. 

In the new algorithm, a significant improvement in efficiency is gained since the production 
coelficients are now obtained through a simple log-linear interpolation and the exponential attenuation 
terms are computed directly from the attenuation cross sections. Once the main loop is entered, only the 

interpolated values, ^.(®,r + v.h ) , must be computed at each iteration. No integral fluences need to be 

computed, and no interpolations over the integral fluences need to be calculated. The most computationally 
expensive portion of the calculation is now the interpolation and has therefore been given a dashed outline 
in Figure 3. The interpolation routine [Wilson et al. 1995] used in previous versions of HZETRN and 
BRYNTRN was reviewed and several modifications were made to improve efficiency and robustness. The 
updated routine still performs log-log cubic interpolation with Lagrange polynomials but is more efficient 
and has improved logic for interpolation around certain discontinuities and extrapolations past r max . A 
comparison of the routines will be given later. It should also be noted that in the new algorithm, all of the 
computations are carried out in double precision; the need for double precision will be shown later. 

The updated algorithm is much simpler and faster than its predecessor. An estimate of the 
improvements can be obtained by simply counting the number of interpolations in a single iteration of the 
light particle transport equation. For this exercise, recall that there are six particles to transport and assume 
there are N grid points in the r-grid. The results are summarized in Table 1. It is clear the new algorithm 
requires several orders of magnitude fewer interpolations than the old algorithm. 
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Figure 3. Process description of the new light particle transport algorithm. 


Table 1. Number of interpolations per step in the old and new light particle transport algorithms. 


Grid-size 

(AO 

Interpolations 

Old 

New 

100 

364,800 

600 

150 

817,200 

900 

300 

3,254,400 

1800 

500 

9,024,000 

3000 


7. Heavy Ion Marching Equation Numerical Methods 

The heavy ion marching equation was given previously as 

ipjix + h,r) = 

+ E— ( r + v i h / 2 )M x ’ r + ( u j + v k) h / 2 ) ( 95 > 

k>j ^ k 

g— °j(r +Vjh/2)h _ e -^ h (T+( v j+ v h)hl*)h 

X a k (r + (v j + v k )h / 2) - a-(r + Vj h / 2) ' 
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Inspection of equation (95) reveals that the only numerical technique required is interpolation. The 
interpolation routine used in the heavy ion marching algorithm [Wilson et al. 1995] is the same routine 
used in the light particle marching algorithm. Thus, to improve efficiency and robustness, this routine has 
been replaced with the updated version mentioned in the previous section. A comparison of the routines 
will be given shortly. 

8. Coupling Light Particle and Heavy Ion Solutions for GCR 


The GCR environment is composed of energetic protons, alpha particles, and heavier ions with 
Z > 2 . In order to simultaneously transport all of these particles using the methods outlined above, the 
light particle and heavy ion marching algorithms must be coupled. This coupling has never been explicitly 
documented, but has existed in HZETRN for some time. The complete marching algorithm for GCR 
environments can be succinctly written as 

ipfa + A,r.) = + "ft) 

+ H ( J - j)E E a imM X > r m) 

k 771=1 

^ v. . (96) 

+ E ~ a jk ( r + v j h / 2 )M x ’ r + ( v i + v k) h / 2 ) 

k>J* V k 

e -<T } (r+v } h/2)h _ e -o t (r+(L' j +v l .)h/2)h 
a k( r + ( v j + v k) h / 2 ) “ a j( r + v j h / 2 ) 

The index J* refers to the heaviest light ion ( 4 He), the Heaviside function is 


H(J - j ) 


0 , j > J 

1 , i < J* 


and the lower limit of the last summation term is given by 


t* = . j ’ j> J . 

J i 3 — J 


(97) 


(98) 


For completeness, the same marching equation is given in terms of the old marching algorithms as 

ipj(x + h,r) = + Vjh) 

+h(j - i)E^E e '" i(ri)V2 " ff * (rJV2i ^( r i^ m + / 2;^) 

k ^k m—i 

x [ * j & r m + 0 v i + v k ) h / 2) - *j (*. r m+1 + (Vj + v k )h / 2) ] (99) 

+ E ~ a jk ( r + v j h / 2 )M x > r + ( u j + v k) h / 2) 

k>J* 

^-o^r+v-hl^h _ e ~(r k {r+(v j +is k )h/2)h 

X a k (r + (Vj + v k )h / 2) - a^r + / 2) ‘ 

The only difference between equation (96) and (99) is the form of the light particle production term. The 
physical meaning of the coupling is the same in either equation. 

The simplest way to explain and understand the GCR transport equation is by considering two 
separate cases. First, suppose we are interested in the transport of a single heavy ion with particle subscript 
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p. Since p refers to a heavy ion, it must be true that p > J , so that H(J — p) = 0 and J* = p . 
Therefore, the light particle production term is neglected and only the heavy ion transport equation remains 

+ h,r i ) = e~ a ^ h 'tp p {x,r i + v p h) 

+Yj— a P k ( r + v P h / 2 )M x ’ r + (y P + v k) h / 2 ) ( 10 °) 

k>p ^ k 

e -<r p (r+v p h/2)h _ e ~(T k {r+{v p +v k )hl2)h 
a k( r + K + v k) h / 2) - a p( r + v p h / 2) ' 

There is nothing surprising in this result; it simply reiterates that heavy ions are transported with the heavy 
ion transport equation. 

Now, suppose that we are interested in the transport of a single light particle with particle 
subscript p. Since p refers to a light particle, it must be true that p < J , so that H(J — p) = 1 and 
J* = J . In this case, the marching algorithm becomes 

+ h,r]) = e~ ar(r>)h tp p {x,r i + u p h) 

+E E a im M x > r J 
k 771=1 

^1/ . (101) 

+E — a P k ( r + u P h / 2 )M x > r + ( u P + v k) h / 2 ) 

k>J V k 

e ~cr p (r+v p h/2)h _ e -<T k (r+(v p +v k )h/2)h 

X(7 k( r + ( V P + v k) h / 2) - + v p h / 2) ’ 

which is identical to the light particle marching algorithm, except that the additional summation term has 
been added. The second summation term is taken over all heavy ions and represents the source of light 
particle projectile fragments produced by collisions between heavy ions and the target media. This second 
summation term is the coupling mechanism between the light particle and heavy ion marching equations. 

9. Round-Off Error, Coding Mistakes, and Interpolation 

In this section, the impact of single precision round-off error, coding errors and interpolation on 
exposure quantities and overall code stability is examined. First, the errors caused by single precision 
calculations in the light particle cross sections are analyzed. These errors are significant enough to cause 
instabilities in the transport algorithms that are problematic for material thicknesses > 50 g/cm 2 . Various 
coding errors related to light particle production cross sections that have a non-negligible impact on fluence 
spectra and exposure quantities are also discussed. Finally, the new interpolation routine mentioned in 
sections 6 and 7 is compared to the previous routine [Wilson et al. 1995], and the improved efficiency and 
robustness of the updated algorithm is shown. For the convergence tests (discussed later), all the errors, 
mistakes, and interpolation routines have been corrected in the original and updated algorithms so that any 
differences can be attributed only to the numerical algorithms. 

Subtractive cancellation can occur in computational algorithms when the difference is computed 
between two numbers that are nearly equal in a given precision. This problem was occurring for some of 
the light particle cross section calculations. Recall that the integrated light particle production cross sections 
were computed as a difference of two cumulative production spectra according to 


F^r^h) = “ l a jk (E, E ') dE - a jk (E,E')dE 


e(r) 


( 102 ) 
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For large fragment energies, where e(r + VjK) ~ s(r), the two integrals in equation (102) become nearly 

equal. Thus, in single precision, the difference will lose numerical precision and may be evaluated as zero. 
An example of this is shown in Figure 4, where the integrated light particle production cross section, 

(r, r h ) , for 1 AGeV 4 He ions producing 3 H and 3 He in aluminum [Cucinotta et al. 1 993] are evaluated 

in single and double precision. For fragment energies larger than ~50 AMeV, the single precision results 
are zero and are therefore not visible on the plot. Similar results were also found for other reactions and 
projectile energies. The impact of this round-off error on particle fluence spectra can be seen in Figure 5, 
where the secondary 3 H and 3 He fluences at 100 g/cm 2 in aluminum exposed to the 4 He component of the 
1977 solar minimum GCR spectrum [O'Neill and Badhwar 2006] are shown. The instabilities in the single 
precision results above 10 AGeV are clear. As these results are propagated to larger depths, the instabilities 
grow in magnitude, reach lower energies, and eventually cause algorithm failure. 

Various coding errors were also discovered in the algorithms that generate the light particle 
production cross sections. To understand these coding errors, it is necessary to first review how differential 
production cross sections are represented in HZETRN. From Wilson et al. [1991], the light particle 
differential cross sections can be expressed as 

a jk (E, E ') = m jk (E ')* k (E % k (E, E ') , (103) 

where a k (E) is the macroscopic cross section for the projectile, m Jk (E) is the mean multiplicity of type j 
particles produced by type k particles with energy E , and f jk (E 9 E) is the energy spectrum of the reaction. 
The spectral term is normalized so that 


f 0 °°f jk (E,E')dE= 1. 


(104) 



Figure 4. Integrated production cross section for a 1 AGeV 4 He ion producing 3 He and 3 H in aluminum. 
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Figure 5. Secondary 3 H and 3 He fluences computed using single and double precision cross sections at 100 
g/cm 2 in aluminum exposed to the 4 He component of the 1977 solar minimum GCR spectrum. 


For each of the reactions, there is a maximum fragment energy, E m f 9 above which the spectral term 
is identically zero. For example, for nucleon production from nucleon-nucleus interactions, the fragment 
energy is required to be less than the projectile energy; thus, E mf = E [Wilson et al. 1991]. These 
constraints were imposed on the integrated light particle production cross sections in equation (102) by 

setting F* (r, r h) = 0 if either e(r + Vjh) or e(r) were greater than E mf . This constraint correctly accounts 

for the case when both energies are greater than E mf9 but it does not account for the case when e(r + Vjh) > 

E mf and s(r) < E mf . One would expect this error to cause a small under-prediction of secondary particles by 
effectively truncating the integration limits in equation (102). The necessary corrections have been made so 
that the integrated light particle production cross sections are now properly evaluated up to their limits for 
all energies and step-sizes. The effect of this algorithmic error is shown in Figure 6 where we have 
computed the proton fluences at 100 g/cm 2 in aluminum exposed to the August 1972 King SPE [King 
1974]. As expected, the results with the maximum fragment energy coding error, "E mf Error," under-predict 
the corrected results, "E mf Corrected.” The numerical errors shown in Figure 6 are on the order of 25% 
below 1 MeV. 

Another coding error was found in which a modified proton elastic cross section was incorrectly 
added into a non-elastic reaction channel. For the reaction, n + T — ► p + X, where n is a neutron, T is the 
initial target state, p is a proton, and X is the final target state, the production cross section was evaluated as 

a pn {E,E') = [m pn (E')o n (E')f pn (E,E')\ + a*J(E')f£(E,E') , (105) 

where the superscript "el" refers to elastic. The bracketed term contains the correct quantities for the 
reaction, and the non-bracketed term is the modified proton elastic cross section that should not be 
included. Notice that the additional term uses the neutron elastic cross section and the proton elastic 
spectrum. Thus, not only is the term used incorrectly, it is being evaluated improperly as well. Surprisingly, 
this error, coupled to the E mf error mentioned above, can produce results very near the corrected code. In 
Figures 7 and 8, the isolated and coupled effects of both errors are shown. In the Figures, "E mf Error" refers 
to the maximum fragment energy error mentioned previously, "p el Error" refers to the proton elastic error in 
equation (105), and the "Corrected" results were obtained with both errors removed. The coding errors have 
very different impacts in aluminum and water. In the aluminum target, the proton elastic error has almost 
no effect when the maximum fragment energy error is present. However, in the water target, the proton 
elastic error has a large effect whether the maximum fragment energy error is present or not. 
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Figure 6. Proton fluences computed with and without the maximum fragment energy coding error at 100 
g/cm 2 in aluminum exposed to the August 1 972 King SPE. 


A final coding error was found in the evaluation of direct knockout target fragments with mass 2 < 
A < 4. The spectral term for this reaction was parameterized as [Cucinotta et al. 1996] 


f jk ( E , E ') 


Aj exp (~EA j / w 0 ) 
w o I 1 — ex P(— ^'/ «» 0 )] ’ 


(106) 


where Aj is the fragment mass and w 0 is a spectral width dependent on the projectile energy given by 


w 0 (E') 


40 + E'/5 , if E' < 800 
200 , if E' > 800 


(107) 


After reviewing the algorithm where this parameterization is found, it was determined that the fragment 
energy, E , instead of the projectile energy, E, was being incorrectly used in the evaluation of w 0 . The effect 
of this spectral width coding error on secondary particle production is shown in Figure 9 where we have 
shown the 4 He and secondary 2 H flux spectra in 100 g/cm 2 of aluminum exposed to the 1977 solar 
minimum GCR environment. The dashed curves marked with "w 0 Error” are the results generated with the 
fragment energy in place of the projectile energy in the evaluation of w 0 in equations (106) and (107). The 
impact is significant between 10 AMeV and 1000 AMeV. Similar results were also found in the secondary 
3 H and 3 He spectra as well. 
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Figure 7. Proton fluences computed with and without the maximum fragment energy and proton elastic 
coding errors at 100 g/cm 2 in aluminum exposed to the August 1972 King SPE. 



Figure 8. Proton fluences computed with and without the maximum fragment energy and proton elastic 
coding errors at 100 g/cm 2 in water exposed to the August 1972 King SPE. 
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Figure 9. 4 He and secondary 2 H flux spectra computed with and without the spectral width error in 100 
g/cm 2 of aluminum exposed to the 1 977 solar minimum GCR environment. 


Lastly, the interpolation routine used in previous versions of HZETRN is compared to a new 
interpolation routine that is faster and more robust. Both routines used cubic Lagrange interpolating 
polynomials. A cubic Lagrange interpolating polynomial requires four data points (jq, yO, (jq+i, y i+ i), (jq+ 2 , 
y i+2 ), and (jq+ 3 , ^+ 3 ). A common choice is to choose the four abscissa points symmetrically about the 
interpolating point so that jq+i < x < jq+ 2 . However, if the interpolating point is near the boundary of the 
domain, one may also choose jq < x < jq+i or jq +2 < x < jq +3 . The log-log cubic Lagrange interpolating 
polynomial is obtained by taking the natural log of the four data points, (jq, y*), (jq+i, y i+ 1 ), (jq+ 2 , y i+2 ), (jq+ 3 , 
y i+ 3 ), and the interpolation point, x. The interpolated value is then exponentiated to account for the initial 
logarithm. 

The interpolation routine previously used in the transport algorithms [Wilson et al. 1995] 
performed log-log cubic interpolation and extrapolation. The search procedure used to find the abscissa 
points closest to the interpolation point was linear and searched in one direction from the origin (started 
from the first data point and searched forward one point at a time). All logarithms and exponentials were 
computed within the algorithm. No logic was included to ensure that extrapolated results matched the 
trends of the data within the domain (ie. increasing or decreasing). Also, no logic was included to check for 
non-smooth data. 

In the updated algorithm, the search procedure used to find the abscissa points is linear and 
searches in either direction from the previously used interpolation point. Thus, when the routine is called 
repeatedly with an ordered or sequential set of interpolation points, the search algorithm uses previous 
results to find the nearest abscissa points. The new algorithm also avoids the use of logarithms and 
exponentials. Instead, the natural log of the data points is computed once prior to interpolation and used 
repeatedly. The extrapolation procedure has been switched to log-log linear. This will ensure that 
extrapolated results always match the trend of the data near the end of the domain. Logic has also been 
included to check for non-smooth data. In the event that the interpolated result does not lay between its 
nearest two data points, the interpolated value is re-computed using linear interpolation. The benefit of this 
logic for non-smooth data sets will be examined shortly. 
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To show the improved efficiency of the updated algorithm, the August 1972 King SPE is 
transported through 100 g/cm 2 of aluminum using a step-size of h = 0.5 g/cm 2 and neglect all nuclear 
production. The transport solution is given simply as 


ip p (x + h,r) = e <Tp ^ h ip p (x,r + h ) , (108) 

requiring only interpolation to compute ip p (x,r + h ) . It took -4.9 milliseconds to reach 100 g/cm 2 using 

the new interpolation routine, and it took -16.3 milliseconds to reach 100 g/cm 2 using the old interpolation 
routine. The new routine is almost 3.5 times faster than the previous one. 

To show the improved robustness of the updated algorithm, we propagate a modified form of the 
August 1972 King SPE through 2 g/cm 2 of aluminum using a step-size of h = 0.5 g/cm 2 . The modified 
spectrum has been adjusted so that the proton fluence above 50 MeV is 1.0 particles/(cm 2 -MeV). The 
results are shown in Figure 10. Note that the original discontinuity at 50 MeV has been shifted to 
approximately 20 MeV because of atomic interactions in the target. The benefit of the added logic for non- 
smooth data is clear. The old algorithm over-predicts the results below 20 MeV, especially near the 
discontinuity; it also sharply under-predicts the results just past the discontinuity. This oscillatory behavior 
is characteristic of cubic polynomials and can be problematic when such errors are propagated to larger 
depths. Conversely, the new algorithm properly interpolates through a discontinuity and introduces no 
oscillatory behavior. 



Figure 10. Comparison of old and new interpolation routines in 2 g/cm 2 of aluminum exposed to a modified 

form of the August 1972 King SPE. 
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To show the combined effects of the errors identified in this section, we compare dose values 
computed with the original code (single precision, E mf error, p d error, w 0 error, and old interpolation) and 
corrected code at 100 g/cm 2 of aluminum and water exposed to the August 1972 King SPE and 1977 solar 
minimum GCR environments. The data are given in Table 2. First, note that the SPE data in aluminum and 
GCR data in aluminum and water are surprisingly similar despite the individual effects of the errors shown 
above. This would indicate that many of the errors are competing (ie. cancelling effects) and shows how 
difficult these errors can be to find and correct. The large differences in the water target exposed to the SPE 
are caused almost entirely by the p el error. The w 0 error causes the code with errors to under-predict doses 
in aluminum. The results in the water target exposed to the GCR show the combined effects of two errors. 
While the p el error would cause an over-prediction, the w 0 error causes an under-prediction, and the 
combined effect is a result very near the corrected code. It should be noted that for different materials or 
environments, the results are unpredictable. 


Table 2. Dose at 100 g/cm 2 in aluminum and water exposed to August 1972 King SPE and 1977 solar 
minimum GCR. The units for the SPE results are cGy, and the units for the GCR results are cGy/day. 



SPE 

GCR 



Aluminum 

Water 

Aluminum 

Water 

Corrected code 
Code with Errors 

3.28 x 10* 
3.21 x 10 1 

6.18 x 10' 2 
1.27 x 1 0' 1 

4.72 x 1 0' 2 
4.22 x 1 O' 2 

3.78 x 1 0' 2 
3.77 x 1 O' 2 


10. Convergence Study for Light Particle Transport in SPE Environments 

In this section, the total discretization error associated with the old and new light particle transport 
algorithms in HZETRN is examined by conducting a detailed convergence analysis in both step-size ( h ) 
and energy grid-size (N). For all comparisons, the round-off errors, coding mistakes, and interpolation 
routines mentioned in the previous section have been fixed in the original algorithms. This will allow a 
direct comparison between the original and updated numerical algorithms. The convergence analysis is 
completed by transporting the August 1972 King SPE spectrum through 100 g/cm 2 of aluminum and water. 
Six different energy grids were used with minimum and maximum energy values of E^n = 0.01 AMeV and 
Ermx = 2500 AMeV. The grids were equally log spaced in energy and contained 100, 149, 223, 335, 502, 
and 753 points; these grids will be referred to as £-100, £-149, £-223, £-335, £-502, and £-753. Each grid 
represents a refinement of 1 .5 in the grid spacing parameter, 

A = lQ g(£,„J - lQ g(^niJ a09) 

E N - 1 


Thus, for N= 100, A E = 0.054525, and for N = 149, A E = 0.036473. Along with these different energy grids, 
the following step-sizes (in g/cm 2 ) were used to propagate the boundary condition into the slabs: h = 2' 1 , 
2' 2 , 2' 3 ,...,2' n . Each step-size is a refinement of two in h. These eleven step-sizes, six energy grids, two 
materials, and two transport algorithms, were used to obtain fluence, dose, and dose equivalent values at 
various depths in the target media. These data are used to show that both algorithms converge as a function 
of step-size and energy grid-size for SPE boundary conditions. It is also determined that the new algorithm 
is more accurate than its predecessor. Finally, total discretization error estimates are given for the 
discretization parameters (h = 0.5 g/cm 2 with 100 energy grid points) commonly used in HZETRN. The 
errors are expressed as percent difference from a highly accurate, or converged, numerical solution. 

In Figures 11 and 12, dose equivalent as a function of depth in aluminum and water is given for 
both of the algorithms and all of the discretization parameters. The plot legends have been removed 
because there are 132 different curves (6 energy grids, 11 step-sizes, 2 algorithms), many of which are 
overlapping. The spread in the data at 100 g/cm 2 in aluminum is -3 cSv (-39%), and the spread in the data 
at 100 g/cm 2 in water is -0.2 cSv (-26%). It is clear in these plots that both methods, with the various 
discretization parameters, produce similar results and that the errors appear to be slowly increasing as a 
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function of depth in the material. Slightly larger errors are also found in aluminum. These errors will be 
discussed in detail later. 

In Figures 13 and 14, the dose equivalent at 100 g/cm 2 in aluminum and water is plotted as a 
function of the energy grid-size, N 9 for three step-sizes. The three step-sizes were chosen to reduce data 
overlap and improve plot clarity. In Figures 15 and 16, the dose equivalent at 100 g/cm 2 in aluminum and 
water is plotted as a function of the step-size, h 9 for all of the energy grids. Figures 13-16 show that both of 
the algorithms reach an asymptotic solution for step-sizes less than 0.01 g/cm 2 and energy grids with 
greater than 300 points. The use of step-sizes larger than 0.01 g/cm 2 results in a systematic under-prediction 
of dose equivalent regardless of the method or energy grid used. It is also clear that if larger energy grids 
and smaller step-sizes were considered, the methods would generate almost identical solutions. The percent 
difference between the dose equivalent values generated by the old and new algorithms with the finest 
discretization parameters (h = 2' 11 , N = 753) is 0.4% at 100 g/cm 2 in aluminum and 0.2% at 100 g/cm 2 in 
water. Thus, it is concluded that both methods converge to the same value as a function of step-size and 
energy grids size to the same solution for this SPE in aluminum and water. 

Figures 15-16 indicate that there is little difference in spatial discretization error between the old 
and new algorithms. This is expected since the only difference between the algorithms is in the calculation 
of the source integral over energy. The new method appears to have a smaller energy discretization error in 
aluminum. This can be quantified by computing the percent difference between the £"-100 and £'-753 dose 
equivalent results from both algorithms with a step-size of h = 2~ u g/cm 2 at 100 g/cm 2 . The percent 
difference between the £"-100 and £- 753 results for the old algorithm is 10.7% in aluminum and 2.4% in 
water, while the percent difference between the £-100 and £-753 results for the new algorithm is 2.5% in 
aluminum and 2.6% in water. The improvement in aluminum was achieved primarily because the new 
algorithm more accurately handles the sharply peaked production cross sections in aluminum, especially 
the neutron elastic cross sections that are important at large depths. The negligible accuracy loss in water 
shows that the methods are comparable when the production cross sections are relatively smooth, as is the 
case when hydrogen is present. 



Figure 11. Dose equivalent versus depth in aluminum exposed to the August 1972 King SPE computed 
with both transport algorithms and all discretization parameters. 


29 




Figure 12. Dose equivalent versus depth in water exposed to the August 1972 King SPE computed with 
both transport algorithms and all discretization parameters. 



Number of energy grid points N 

Figure 13. Dose equivalent versus energy grid-size at 100 g/cm 2 in aluminum exposed to the August 1972 
King SPE. Both transport algorithms were used with three different step-sizes. 
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Number of energy grid points N 


Figure 14. Dose equivalent versus energy grid-size at 100 g/cm 2 in water exposed to the August 1972 King 
SPE. Both transport algorithms were used with three different step-sizes. 



Figure 15. Dose equivalent versus step-size at 100 g/cm 2 in aluminum exposed to the August 1972 King 
SPE. Both transport algorithms were used with six different energy grids. 
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Figure 16. Dose equivalent versus step-size at 100 g/cm 2 in water exposed to the August 1972 King SPE. 
Both transport algorithms were used with six different energy grids. 


While one would like to use the fully converged results in all future applications, the fine 
discretization parameters result in a very inefficient algorithm that would be impractical in most 
applications. Thus, a smaller grid and larger step-size is chosen, and the results generated with these coarse 
discretization parameters must be compared to the converged results to quantify the discretization error. As 
stated previously, HZETRN is commonly run with 100 energy grid points and a step-size of 0.5 g/cm 2 . 
Therefore, h = 0.5 g/cm 2 and N = 100 are chosen as the coarse discretization parameters. Define the 
converged solution as the results obtained with the new method using 753 energy grid points and a step- 
size of 2' 11 g/cm 2 . This approximation is justified in Figures 13-16 where it appears that both methods are 
asymptotically approaching a common solution. It is also important to note that the difference between the 
methods with the fine discretization parameters is very small. The discretization error is expressed as the 
percent difference between the results obtained with the coarse parameters and the converged solution. The 
results are given in Figures 17 and 18. The new algorithm has a lower discretization error than the old 
algorithm out to 100 g/cm 2 . Figures 17 and 18 also show that characteristically different errors are 
generated in aluminum and water targets. The data represented in Figures 17 and 18 are also given in Table 
3 so that the exact numbers are clearly documented. Table 4 gives the error estimates for dose values. 

Notice that in Figures 17 and 18, there is a slight bend in the error curves near 40 g/cm 2 . This bend 
is the approximate depth at which the contribution to dose equivalent by the primary protons is overtaken 
by the contribution from secondary neutrons and charged target fragments, as shown in Figures 19 and 20. 
In aluminum, the error curve continues to increase past the bend, while in water, the error curve decreases 
past the bend. This is due to the number and energy of secondary neutrons and charged target fragments 
produced in each material. In aluminum, there are a moderate number of neutrons produced as a result of 
nuclear collisions between the primary protons and target media, and the many elastic collisions which 
dominate neutron transport result in only small energy transfers (compared to neutron collisions with 
hydrogen atoms in water). 
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Figure 17. Discretization error for dose equivalent as a function of depth for the light particle transport 
algorithms in aluminum exposed to the August 1972 King SPE. 



Figure 18. Discretization error for dose equivalent as a function of depth for the light particle transport 
algorithms in water exposed to the August 1972 King SPE. 
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Table 3. Discretization error (%) for dose equivalent computed with the old and new light particle transport 
algorithms at various depths in aluminum and water exposed to the August 1972 King SPE. Significant 
figures shown to clarify magnitudes. 


Depth (g/cm 2 ) 

Old 

Aluminum 

New 

Old 

Water 

New 

1.0 

0.48 

0.49 

0.57 

0.61 

5.0 

0.08 

0.10 

0.95 

0.97 

10.0 

0.71 

0.63 

2.23 

2.14 

25.0 

10.10 

8.94 

15.24 

14.01 

50.0 

23.40 

18.77 

23.83 

20.73 

75.0 

27.42 

20.18 

22.43 

18.65 

100.0 

31.48 

21.82 

20.80 

16.67 


Table 4. Discretization error (%) for dose computed with the old and new light particle transport algorithms 
at various depths in aluminum and water exposed to the August 1972 King SPE. Significant figures shown 
to clarify magnitudes. 


Depth (g/cm 2 ) 

Old 

Aluminum 

New 

Old 

Water 

New 

1.0 

0.09 

0.79 

0.31 

0.37 

5.0 

0.03 

0.04 

0.38 

0.44 

10.0 

0.11 

0.11 

0.49 

0.51 

25.0 

1.63 

1.45 

2.93 

2.71 

50.0 

17.49 

14.56 

15.43 

13.65 

75.0 

26.53 

20.33 

16.20 

13.53 

100.0 

30.54 

22.04 

14.51 

11.46 


Thus, the secondary neutrons are left with sufficient energy to continue producing charged target fragments 
well past 40 g/cm 2 . Conversely, in water, there are far fewer neutrons produced, and the elastic collisions 
result in much larger energy transfers that render many of the neutrons incapable of producing further 
target fragments. It is important to note that the error curves in Figure 17 do not grow without bound; they 
reach a maximum of about 75% near 725 g/cm 2 and then decline rapidly as the remaining neutrons have 
insufficient energy to produce any further target fragments. Error estimates for material thicknesses greater 
than 100 g/cm 2 are computationally expensive due to the long run-times associated with fine discretization 
parameters and large depths. This topic will be investigated in future work. 

As stated previously, it is clear in Figures 15 and 16 that both methods, for a given energy grid, 
systematically under-predict dose equivalent values if step-sizes larger than 0.01 g/cm 2 are used. It is now 
shown that these step-sizes, including the coarse discretization parameter, h = 0.5 g/cm 2 , commonly used in 
HZETRN, result in a systematic under-prediction of secondary target fragments, and that this under- 
prediction is the dominant source of error shown in Figures 17 and 18. 
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Figure 19. Dose equivalent by particle type using the converged results at 100 g/cm 2 in aluminum exposed 

to the August 1972 King SPE. 



Figure 20. Dose equivalent by particle type using the converged results at 100 g/cm 2 in water exposed to 

the August 1972 King SPE. 
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Recall from Section 5 that two convergence criteria were identified in developing numerical 
procedures for the light particle transport equation, namely that 

(HO) 

a j ( r ) 

and 

h < r-. (Ill) 

The first criterion (step-size much smaller than the nuclear mean free path) is trivially satisfied since 
nuclear mean free path lengths are typically on the order of many cm, but the second criterion (step-size 
much smaller than the range), which has been derived in the present work, is not satisfied for low energy 
particles unless h is taken several orders of magnitude below 0.5 g/cm 2 . 

To show the error induced by choosing discretization parameters that violate equation (111), a 
portion of the convergence test outlined previously is re-run with the minimum energy increased to 50 
AMeV. This energy was chosen because 50 MeV protons have a range ~2.9 g/cm 2 , which should allow 
step-sizes of 0.5 g/cm 2 to be used with only minimal error. Since the systematic under-prediction has been 
shown to be independent of the numerical algorithm, only the new light particle transport equation was 
used. Also, since the error was most prevalent in aluminum, the water target was not considered. For this 
comparison, define the modified converged solution as the results generated by the new method with a 
step-size of 2' 11 g/cm 2 and the 753 point modified energy grid (iT m j n = 50 AMeV instead of 0.01 AMeV). 
The results generated by the new method using a step-size of 0.5 g/cm 2 and the modified 100 point energy 
grid will be compared to the modified converged results. The error comparisons are given in Figure 21 
along with the previous error curve shown in Figure 17. It is clear that neglecting the target fragments with 
energy less than 50 AMeV significantly reduces the discretization error in the light particle transport 
algorithm. The new error curve has a maximum of -0.07% near 30 g/cm 2 of aluminum. 



Figure 21. Discretization error (%) caused by target fragments with energy less than 50 AMeV in 
aluminum exposed to the August 1972 King SPE. 
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11. Convergence Study for HZETRN in GCR Environments 


For GCR environments, the primary spectrum is composed of protons, 4 He, and 26 heavy ions 
with energies extending up to 50 AGeV. In HZETRN, the heavy ion transport algorithm is used to transport 
the 26 heavy ions along with their secondary fragments through the target media, and the light particle 
transport algorithm is used to transport the primary protons and 4 He along with the light projectile and 
target fragments. The light projectile fragments produced from the heavy ions are also included in the light 
particle transport algorithm. A convergence analysis has never been done, to our knowledge, for the heavy 
ion or light particle transport algorithms in GCR environments. The convergence analysis described in the 
previous section is repeated for the heavy ion and light particle transport algorithms for GCR boundary 
conditions. 

In this section, the total discretization error associated with the heavy ion transport algorithm and 
the old and new light particle transport algorithms in HZETRN is examined by conducting a detailed 
convergence analysis in both step-size, h 9 and energy grid-size, N. For all comparisons, the round-off 
errors, coding mistakes, and interpolation routines mentioned in the previous section have been fixed in the 
original algorithms. This will allow a direct comparison between the original and updated numerical 
algorithms. The convergence analysis is identical to the one in Section 10 except that 1977 solar minimum 
GCR environment is used, and the maximum energy is increased to E max = 50 AGeV. The same step-sizes 
and energy grid-sizes are used. It is first shown that the heavy ion transport algorithm converges as a 
function of step-size and energy grid-size. The discretization error associated with the heavy ion transport 
algorithm for the discretization parameters (h = 0.5 g/cm 2 with 100 energy grid points) commonly used in 
HZETRN is also given. Next, the two light particle transport algorithms are compared. The new algorithm 
is shown to converge as a lunction of step-size and energy grid-size. It is also determined that the old 
algorithm converges as a function of step-size, but that larger energy grids are required to determine if the 
old algorithm converges as a function of energy grid-size. Finally, total discretization error estimates are 
given for discretization parameters commonly used in HZETRN. The errors are expressed as percent 
difference from a converged numerical solution. 

In Figures 22 and 23, dose equivalent as a function of depth in aluminum and water is given for 
the heavy ion transport algorithm coupled to both of the light particle transport algorithms and all of the 
discretization parameters. The plot legends have been removed because there are 132 different curves, 
many of which are overlapping. The spread in the data at 100 g/cm 2 in aluminum is -0.024 cSv/day 
(-17%), and the spread in the data at 100 g/cm 2 in water is -0.013 cSv/day (-15%). As before, both 
methods with the various discretization parameters produce similar results, and the errors appear to be 
slowly increasing as a function of depth in the material. Slightly larger errors are once again found in 
aluminum. In Figures 24 and 25, the contribution to dose equivalent from heavy ions at 100 g/cm 2 in 
aluminum and water is plotted as a function of the energy grid-size, A, for five step-sizes. The step-sizes 
were chosen to reduce data overlap and improve plot clarity. In Figures 26 and 27, the contribution to dose 
equivalent from heavy ions at 100 g/cm 2 in aluminum and water is plotted as a function of the step-size, h 9 
for all of the energy grids. Figures 24-27 clearly show that the heavy ion marching algorithm reaches an 
asymptotic solution for step-sizes less than 0.01 g/cm 2 and energy grids with greater than 250 points. It is 
concluded that the heavy ion transport algorithm converges as a function of step-size and energy grid-size. 

In order to quantify the discretization error associated with the heavy ion transport algorithm, 
consider the contribution to dose equivalent from heavy ions only and define the converged heavy ion 
solution as the results obtained using 753 energy grid points and a step-size of 2" 11 g/cm 2 . The step-size, h = 
0.5 g/cm 2 , and energy grid-size, N = 100, are once again chosen as the coarse discretization parameters. 
The discretization error is expressed as the percent difference between the results generated with the coarse 
parameters and the converged solution. To be clear, these discretization error estimates were obtained using 
only the heavy ion contribution to dose equivalent. Thus, the error caused by the light particle transport 
algorithm is not included here. The results are shown in Figure 28 and 29 for the aluminum and water 
targets. The discretization errors are bounded over the depths shown by 0.3% in aluminum and 1.1% in 
water. Future work will focus on quantifying the discretization error depths greater than 100 g/cm 2 . The 
data represented in Figures 17 and 18 are also given in Table 5 so that the exact numbers are clearly 
documented. Table 6 gives the error estimates for dose values. 
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Figure 22. Dose equivalent versus depth in aluminum exposed to the 1977 solar minimum GCR computed 
with both transport algorithms and all discretization parameters. 



Figure 23. Dose equivalent versus depth in water exposed to the 1977 solar minimum GCR computed with 
both transport algorithms and all discretization parameters. 
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Figure 24. Dose equivalent from heavy ions versus energy grid-size at 100 g/cm 2 in aluminum exposed to 
the 1977 solar minimum GCR environment. Five different steps-sizes were used. 



Number of energy grid points N 


Figure 25. Dose equivalent from heavy ions versus energy grid-size at 100 g/cm 2 in water exposed to the 
1977 solar minimum GCR environment. Five different steps-sizes were used. 
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Figure 26. Dose equivalent from heavy ions versus step-size at 100 g/cm 2 in aluminum exposed to the 1977 
solar minimum GCR environment. All of the different energy grids were used. 



Figure 27. Dose equivalent from heavy ions versus step-size at 100 g/cm 2 in water exposed to the 1977 
solar minimum GCR environment. All of the different energy grids were used. 
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Figure 28. Discretization error for dose equivalent as a function of depth for the heavy ion transport 
algorithm in aluminum exposed to the 1977 solar minimum GCR environment White space included so 

vertical axis matches Figure 29. 



Figure 29. Discretization error for dose equivalent as a function of depth for the heavy ion transport 
algorithm in aluminum exposed to the 1977 solar minimum GCR environment. 
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Table 5. Discretization error (%) for the heavy ion contribution to dose equivalent computed with the heavy 
ion transport algorithm at various depths in aluminum and water exposed to the 1977 solar minimum GCR. 
Significant figures shown to clarify magnitudes. 


Depth (g/cm 2 ) 

Aluminum 

Water 


0.0026 

0.037 


0.024 

0.18 


0.045 

0.32 











0.25 

1.01 


Table 6. Discretization error (%) for the heavy ion contribution to dose computed with the heavy ion 
transport algorithm at various depths in aluminum and water exposed to the 1977 solar minimum GCR. 
Significant figures shown to clarify magnitudes. 


Depth (g/cm 2 ) 

Aluminum 

Water 

1.0 

0.000327 

0.020 

5.0 

0.017 

0.13 

10.0 

0.035 

0.23 

25.0 

0.070 

0.43 

50.0 

0.12 

0.61 

75.0 

0.16 

0.73 

100.0 

0.20 

0.82 


In Figures 30 and 31, dose equivalent at 100 g/cm 2 in aluminum and water is plotted as a function 
of the energy grid-size, N, for three step-sizes. The step-sizes were chosen to reduce data overlap and 
improve plot clarity. In Figures 32 and 33, dose equivalent at 100 g/cm 2 in aluminum and water is plotted 
as a function of the step-size, h , for all of the energy grids. Figures 30-33 clearly show that the new light 
particle transport algorithm coupled to the heavy ion transport algorithm reaches an asymptotic solution for 
step-sizes less than 0.01 g/cm 2 and energy grids with greater than 300 points. It is concluded that the new 
light particle transport algorithm converges as a function of step-size and energy grid-size. Conversely, 
while the old algorithm appears to converge as a function of step-size, it has not converged as a function of 
energy grid-size if 753 grid points are used. This does not mean that the old algorithm will not converge. 
All that can be said is that an asymptotic solution is not achieved by the old algorithm if as many as 753 
energy grid points are used. This would indicate that the new method is more accurate and converges faster 
as a function of energy grid-size. 

In order to quantify the discretization error associated with the algorithms, define the converged 
solution as the results obtained with the new light particle transport algorithm (coupled to the heavy ion 
algorithm) using 753 energy grid points and a step-size of 2' 11 g/cm 2 . The step-size, h = 0.5 g/cm 2 , and 
energy grid-size, N= 100, are once again chosen as the coarse discretization parameters. The discretization 
error is expressed as the percent difference between the results generated with the coarse parameters and 
the converged solution. These error estimates were obtained using total exposure quantities (heavy ions and 
light particles included), thus the discretization error from the heavy ion transport algorithm, though small, 
is included. The results are shown in Figure 34 and 35 for the aluminum and water targets. The 
discretization errors are bounded over the depths shown by 6.5% in aluminum and 3.9% in water. The 
errors shown for the old method are misleadingly small. First, recall that the old method has not reached an 
asymptotic solution if 753 energy grid points are used; thus, it converges more slowly than the new 
algorithm as a function of energy grid-size. Second, the errors reported here are dependent on the choice of 
the converged solution. In this case, the energy discretization error and the spatial discretization error in the 
old algorithm are competing errors and therefore produce misleading error estimates when compared to a 
converged solution. 
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Figure 30. Dose equivalent versus energy grid-size at 100 g/cm 2 in aluminum exposed to the 1977 solar 
minimum GCR environment. Both transport algorithms were used with three different step-sizes. 
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Figure 31. Dose equivalent versus energy grid-size at 100 g/cm 2 in aluminum exposed to the 1977 solar 
minimum GCR environment. Both transport algorithms were used with three different step-sizes. 
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Figure 32. Dose equivalent versus step-size at 100 g/cm 2 in aluminum exposed to the 1977 solar minimum 
GCR environment. Both transport algorithms were used with six different energy grids. 



Figure 33. Dose equivalent versus step-size at 100 g/cm 2 in aluminum exposed to the 1977 solar minimum 
GCR environment. Both transport algorithms were used with six different energy grids. 
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Figure 34. Discretization error (%) for dose equivalent as a function of depth for both algorithms in 
aluminum exposed to the 1977 solar minimum GCR environment. 



Figure 35. Discretization error (%) for dose equivalent as a function of depth for both algorithms in water 
exposed to the 1977 solar minimum GCR environment. 
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As before, the errors are slowly increasing as a function of depth in the material, and it is expected that 
these errors will reach a maximum after the majority of the primary ions and secondary projectile 
fragments have come to rest. Future work will focus on quantifying the discretization error depths greater 
than 100 g/cm 2 . The data represented in Figures 34 and 35 is also given in Tables 7. Table 8 gives the error 
estimates for dose. 


Table 7. Total discretization error (%) for dose equivalent computed with both algorithms at various depths 
in aluminum and water exposed to the 1977 solar minimum GCR. Significant figures shown to clarify 
magnitudes. 


Depth 

(g/cm 2 ) 

Aluminum 

Old 

New 

Water 

Old 

New 

1.0 

0.13 

0.12 

0.17 

0.17 

5.0 

0.15 

0.28 

0.31 

0.47 

10.0 

0.14 

0.52 

0.39 

0.84 

25.0 

0.07 

1.46 

0.37 

1.88 

50.0 

0.40 

3.27 

0.28 

3.03 

75.0 

1.37 

4.96 

0.63 

3.62 

100.0 

2.92 

6.47 

1.42 

3.97 


Table 8. Total discretization error (%) for dose computed with both algorithms at various depths in 
aluminum and water exposed to the 1977 solar minimum GCR. Significant figures shown to clarify 
magnitudes. 


Depth 

(g/cm 2 ) 

Aluminum 

Old 

New 

Water 

Old 

New 

1.0 

0.04 

0.07 

0.07 

0.11 

5.0 

0.16 

0.12 

0.10 

0.24 

10.0 

0.48 

0.18 

0.41 

0.33 

25.0 

1.53 

0.32 

1.36 

0.46 

50.0 

3.09 

0.56 

2.49 

0.55 

75.0 

4.25 

0.83 

2.80 

0.56 

100.0 

4.91 

1.11 

2.39 

0.53 


12. Run Time Comparisons 

In this section, the run-times of the old and new transport algorithms are compared for three 
applications in which HZETRN is commonly used. All simulations were run on a Dell with 4 quad-core 
AMD Opteron Processor 8356 chips and 131 GB of memory. The discretization parameters used in each 
case where h = 0.5 g/cm 2 and 100 energy grid points. The simplest of the three applications would be to 
use HZETRN to compute flux/fluence, dose, or dose equivalent as a function of depth in a single layer of 
some material. This type of calculation might be done to analyze the radiation protection properties of a 
slab of some material or to generate a simplified interpolation database. In this case, consider separately 
aluminum exposed to the August 1972 King SPE and the 1977 solar minimum GCR spectrum. For the SPE 
environment, results are stored at 21 depths from 0.0 g/cm 2 to 100.0 g/cm 2 . For the GCR environment, 
results are stored at 11 depths from 0.0 g/cm 2 to 100.0 g/cm 2 . The difference in the number of depths has 
been discussed by Anderson et al. [2007] and is related to the rapid decline of the SPE dosimetric quantities 
over the first few g/cm 2 of shielding. The run-times are given in Table 9. The next application is a two- 
layer simulation in which results are stored on a two dimensional grid of points covering every combination 
of first layer thicknesses followed by second layer thicknesses. In this case, consider a first layer of 
aluminum followed by a second layer of water. This simulation is typically used to model a shielding 
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structure (aluminum) surrounding a human phantom (water). The spatial grids mentioned above for the 
SPE and GCR environments are used in this application. The run-times are given in Table 10. The final 
application is a three-layer simulation in which results are stored on a three dimensional grid of points 
covering every combination of first, second, and third layer thicknesses. In this case, consider a first layer 
of aluminum, a second layer of polyethylene, and a third layer of water. This simulation is typically used 
to model a shielding structure (aluminum) with a secondary shield (polyethylene) surrounding a human 
phantom (water). The spatial grids mentioned above for the SPE and GCR environments are used in this 
application. The run-times are given in Table 11. 


Table 9. Single layer run time comparisons for old and new algorithms in SPE and GCR environments. The 
reduction factor is the ratio of the old method run-time to the new method run-time. 



SPE time (sec) 


Old Method 

62.0 


New Method 

0.7 


Reduction Factor 

88.5 

8.5 


Table 10. Two layer run time comparisons for old and new algorithms in SPE and GCR environments. The 
reduction factor is the ratio of the old method run-time to the new method run-time. 



SPE time (sec) 


Old Method 

1272 

926 

New Method 

15 


Reduction Factor 

85 

9 


Table 1 1 . Three layer run time comparisons for old and new algorithms in SPE and GCR environments. 
The reduction factor is the ratio of the old method run-time to the new method run-time. 



SPE time (min) 

GCR time (min) 

Old Method 

440 

173 

New Method 

5 

21 

Reduction Factor 

88 

8 


The run-times in Tables 9-11 show that the new method is much faster than the old method in both 
SPE and GCR environments. The new method is nearly 90 times faster for SPE simulations and nearly 10 
times faster for GCR simulations. It should be noted that run-times in the old method are dependent on the 
atomic complexity of the target material. For targets with many atomic species, the cross sections computed 
within the transport procedure will require more time and therefore slow the transport procedure. However, 
in the new method, cross sections are pre-computed, and therefore, increased atomic complexity does not 
affect the run times. Since the materials considered here had at most two atomic species (water, 
polyethylene), the run times for the old method and reduction factors are both conservative. 

13. Conclusions and Future Work 

In this paper, a detailed derivation of the numerical marching algorithms used for light particles 
and heavy ions in HZETRN was given. Components of the derivation that were previously absent in the 
documentation were discussed in detail. This work serves as a complete set of documentation for the old 
and new numerical methods used in the straight ahead transport procedures in HZETRN. A new numerical 
method was presented for the light particle transport algorithm. The primary differences between the old 
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and new algorithms were summarized in process description charts and tables. In general, the new method 
requires fewer interpolations per step. 

The impact of single precision round-off error and various coding errors on code stability and 
accuracy was discussed. Round-off error was shown to be problematic for some of the light particle 
production cross sections. These errors produce discontinuities in the high energy portion of secondary 
light particle fluence spectra. Three coding errors related to the light particle production cross sections were 
discussed, and the individual and combined effect of the errors on fluence spectra was shown. A new 
interpolation routine used in both the heavy ion and light particle transport algorithm was presented, and 
the improved efficiency and robustness of the routine was demonstrated. 

A convergence study was conducted to determine if the old and new algorithms converge as a 
function of step-size and energy grid-size and to quantify the discretization error associated with both 
algorithms. The old and new algorithms were shown to converge to the same result for the SPE 
environment, and the improved accuracy of the new method was quantified. Total discretization error on 
dose equivalent at 100 g/cm 2 was reduced from 31% to 21% in aluminum and from 20% to 15% in water 
for commonly used energy grids. For the GCR environment, the heavy ion transport algorithm was first 
shown to converge as a function of step-size and energy grid-size, and discretization error estimates were 
bounded over the depths shown by approximately 1% at 100 g/cm 2 in aluminum and water. The new light 
particle transport algorithm was also shown to converge as a function of step-size and energy grid-size. The 
old algorithm was shown to converge as a function of step-size, and it was concluded that energy grids with 
more than 753 points are required to show convergence in the energy domain. This indicates that the new 
method converges faster and is more accurate in GCR environments. Total discretization error on dose 
equivalent at 100 g/cm 2 was 6 % in aluminum and 4% in water for the new method. 

The importance of low energy target fragments on discretization error was clearly shown by 
neglecting all particles with energy less than 50 AMeV in the transport procedure. The error estimates for 
the cut-off energy grid were negligible compared to the full energy grid. This indicates that the dominant 
source of discretization error in HZETRN is caused by low energy target fragmentation. Further, to the 
author's knowledge, only the convergence criterion in equation (110) (step-size smaller than the nuclear 
mean free path length) has been addressed in the literature. The second convergence criterion in equation 
(111) (step-size smaller than the fragment range), found here in the derivation of the numerical procedures, 
has never been explicitly written until now. This work definitively shows the presence of two convergence 
criteria, and that previous convergence studies never addressed the second criterion. 

Run time comparisons between the old and new algorithms showed large improvements for three 
applications in which HZETRN is commonly used. The new algorithm was found to be approximately 90 
times faster for SPE simulations and approximately 10 times faster for GCR simulations. The new 
algorithms will be implemented in all future versions of HZETRN. 

Future work should focus on further reducing discretization error in HZETRN by developing more 
robust methods for handling the transport of low energy target fragments. While physics modeling errors 
associated with the straight ahead approximation may exceed the current numerical error at depths past 100 
g/cm 2 , there is still a need to reduce the numerical error as long as the current model is used in studies with 
such large material depths (atmosphere, lunar surface, martian surface, etc.). 
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16. Appendix A - Inversion of Light Particle Transport Equation 


In this section, the scaled light particle transport equation (equation (23) in Section 3) 

ipj(x,r) = ^2zr s jk (r,r')ip k (x,r')dr' , 


d d , . 

v 4 V cr.fr) 

dx 3 dr 3 


(112) 


is inverted to obtain a Volterra type integral equation. 

For light ions only (neutrons will be handled separately later), define the characteristic variables 




r - i >-x , 


and 


Zj = r + VjX , 


(113) 


(114) 


The differentials with respect to x and r on the left hand side of the light particle transport equation are 
transformed as 


d_ _ _d_^h _d_^_ _ v 

dx dr]- dx d^. dx 1 


_d d 

d Zj dr ljj 


and 


d_ = _d_^Ji , 9 t d 

dr drj. dr d^ dr d^ dr] j 
If we define the following quantities in terms of the characteristic variables 




and 


(115) 


(116) 


(117) 


(118) 


IS ■ noo 

G j (Vj » (j) = J s jk ( r > r ')A {x,r')dr' 


(119) 


the light particle transport equation can be written as 




( 120 ) 


where the division by v . is well defined since is. > 0 for charged particles. Equation (120) can be solved 
using the integrating factor 
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0j(Vji€j) = ~~ P^(» l',Zj)dri', 

3 3 3 2u, J % J J 


where rj Q is a free parameter to be chosen later. Equation (120) is now simplified as 






which can be inverted by integrating from rj Q to rj . . This integration produces 


X.iVi-f,) = --i- P M* . 


2v. J v Q 


This equation must now be transformed back into the original variables (x, r) . The kernel, 
G, (j) f .) , is expressed in terms of the original variables by noting that 


k k 2 




-^f, 00 . r + *,.* + !/• , r + u jX -r, 

i ,r P ’ 1 


Since the fluence at x = 0 is known (see the boundary condition in equation (3)), select the free parameter 
7) q = ^ to obtain 

XjiVo >€j) = = ^(0,Cj) = ^(0,r + i/ y ;c) . (125) 

The integrating factor $j(Vji€j) is expressed in terms of the original variables (x, r) by noting that 
7 y 0 = f has already been chosen, so that 

a / l nj * , . A v , . 


4-M;) = ^rf?*jW>W dl i 


2v. J v 0 

l 


j J 

1 pr—i/jX r + U.X + T ] 1 

= I O' • 

2^ . ^ r+1/ x 3 2 

J 5 ^ 

= ~fo a ^ r + V i t)dt ■ 


Define 


PjM = f’afr + ufldt. 
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The exponential terms appearing in the integrand of equation (123) can also be written in terms of the 
original variables by again noting that r) 0 = f . has been chosen. In this case, 


U j 3 




nr -v.x 
/ a i 

r + v.x + 7 } " 

1 r 7 ?* 

dr)" f <Ti 

r + v.x + rj " 

J r+ t/jX J 

2 

2v . J r+v.x J 
3 3 

2 


nx nx 

= - I <Tj(r + vfidt + l^-rff/r + v.t)dt 




r) *+ v ,x — r 


2v. 


dr}" 


(128) 


Then equation (123) can be written in terms of the original variables ( x , r) as 


tpAx,r) = e ^ r ^ip.(0 ,r + v.x) 


'j k u k 


nr—L/jX noo &j\ 

ri'+vp—r 


r ’ 2i /, 


1 fr+v-x+rj' 6 

J r+v.x *3 

3 2 

( 3 ) 

S jk 


r + v.x + r ] 1 




r + v.x — T } 1 




2 

dr' dr] 1 


(129) 


Let 


rj' = r — v-x + 2v-x 1 , 


(130) 


then equation (129) becomes 


ipjfar) = e W r '%(0 y r + v^) 


+y f f e ^ 3 ^ r,x ^s. k (r + v.x',r')^ k (x — x' y r')dr'dx' 

y JO J r+ VjX' 3 3 


(131) 


At this time, return to equation (23) and seek to obtain a similar expression for neutrons. For 
neutral particles, the continuous slowing down term is identically zero, and so the relevant transport 
equation simplifies to 


'•PnM = 52-=~ n S nk( r ’ r ')M X ’ r ') dr ' ' ( 132 ) 

jfc V k Jr 

The integrating factor defined in equation (127) can be used to invert this equation. Note that for neutrons, 
equation (127) takes the special form 


tA^ e) 


P n (r,x) = xa n (r). 

Equation (132) is expressed in terms of this integrating factor as 


(133) 
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-^\'ip n {x,r)e° n(r) \ = Y^^e xan{r) j°° s nk {ry)^ k {x,r') dr' . 

nr . . 7A r 


Integrate this equation from the boundary at z = 0 to some arbitrary point in the material at z = x , then 
equation (134) becomes 


i>n(x,r) = e * ffn(r V i (0,r) + ^^- f f e°» (r)( * x) s nk {r,r ')ip n (z,r')dr' dz . (135) 

L V, «/ r 


z = x — x\ 


then equation (135) becomes 


ip n (x,r) = e x<Tn{r) rpMr) + f* f°°e (r)a: (y, r - x',r')dr'dz . (137) 

T^ J0 Jr 


Upon further inspection, if j = n in equation (131), then v n = 0, and equation (137) is 

recovered. Therefore, the light particle transport equation can be completely expressed for neutrons and 
light ions in terms of the Volterra type integral equation as 

ipfar) = e“^ (r ’*)^.(0,r + v.x) 

E V • nx poo of (138) 

zr I I e ^ r ’* s ik ( r + v i x \r ] )rp k (x — x\r')dr'dx } . 

^ i/^ Jo J r+ v.x' 3 J 

Equation (138) is the light particle Volterra type integral equation for which numerical procedures will be 
developed. 

17. Appendix B - Inversion of Heavy Ion Transport Equation 

In this section, the scaled heavy ion transport equation (equation (38) in Section 4) 

4~~ v i4~ + a j w (** r ) = a jk ( r )A (*> r ) > ( i39 > 

° x ° r \ k>j v k 

is inverted to obtain a Volterra type integral equation. 

Define the characteristic variables 


V j =r-v j x, 


Zj= r + V j X > 


along with the functions 


XjiVptj) = ipjfrr), 
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VjiVptj) = *j(r)* 


(143) 


and 


VjkiVjitj) = *(*•), (144) 

so that the differential operator on the left hand side of equation (139) transforms as in the previous section, 
and one obtains 

ir-Xjhptj) - = —^-'Z2—v j k(Vk’Qx k (v k ,Q- (145) 

dr )j 2v. 2v jk>j v k 

Equation (145) can be solved using the integrating factor 


= ~Z~ P tjiv'ytjW > 

J J J 2u. J Vo 


(146) 


where r] 0 is a free parameter to be chosen later. Equation (145) is now simplified as 


dr >3 


( ^ )] = ( ^,4 k ) Xk ( Vk ,S k ), (147) 

J 2 ^' k>j^k 


which can be inverted by integrating from rj 0 to r \. . This integration produces 


(%,€,) -^~E-P ^ ) Xk (v'^W . (148) 

2 v jk>j v k Jlk > 


Since the fluence at x = 0 is known (see the boundary condition in equation (3)), we select the free 
parameter rj 0 = f . to obtain 


XjiWo’tj) = 


' C - Vo t + Vo 




= q(0,$) = tfj(0,r + V.x) 


(149) 


The integrating factor $j(Vji€j) is expressed in terms of the original variables (x, r) by noting that 
tJq = f . has already been chosen, so that 
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= irIe^ iT) '^ )dri ' 


1 nr—i/jX 

= J Oi 

2v . ** r+i/.x 3 

3 J \ 

= -fo a t r + v i t)dt 

= -0i(r,x) . 


r + v-x + r ] ' 


F f 


(150) 


The exponential terms appearing in the integrand of equation (148) can also be written in terms of the 
original variables by again noting that rj 0 = £ . has been chosen. In this case, 




1 nr -v-x 

— I <7 

2v ■ J T+ViX 3 
3 J 


2v.Jt 
3 

r + v.x + r ) M 


1 rw' 

drj" f a. 

2u . ” r+VrX 3 
3 3 


r + v.x + rj " 


dr/" 


2i/, 


fo a ^ r + v fi dt + Jv^fzL a j( r + v ] t ) dt 

77 'H- -x — r 1 

^ n 2,. 


(151) 


Then, equation (148) can be written in terms of the original variables as 

ipj(x,r) = e -/ V r -*)^o, r + v.x) 

-j- yp- r^e 

2v j k>J v k Jr+V i x 


- 0 , 


rj'+i/jX—r 


2i/- 


jk 


r + + 7? ' 




r + z^a; — rj r + z^z + r ] ' 


2i/. 


drj 1 . 


(152) 


Define 


r/' = r — z^x + 2^.x' , (153) 

then, equation (152) becomes 

^■(«,r) = e _ ^ (r,ar Vj(0,r + z^z) 

+5^ — f + i/jx'^fy^x — x',r + i/jX^dx' . ^ ^ 

k>j V k 

Equation (154) is the heavy ion Volterra type integral equation for which numerical procedures will be 
developed. 
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